function [Q, M] = getQM(n_seg, n_order, ts)
    Q = [];
    M = [];
    M_k = getM(n_order);
    for k = 1:n_seg
        %#####################################################
        % STEP 2.1 calculate Q_k of the k-th segment 
        T = ts(k);
        % the sequence is p0,p1,...,pn
        Q_k = [0, 0, 0, 0,        0,         0,         0,          0;
               0, 0, 0, 0,        0,         0,         0,          0;
               0, 0, 0, 0,        0,         0,         0,          0;
               0, 0, 0, 0,        0,         0,         0,          0;
               0, 0, 0, 0,  576*T^1,  1440*T^2,  2880*T^3,   5040*T^4;
               0, 0, 0, 0, 1440*T^2,  4800*T^3, 10800*T^4,  20160*T^5;
               0, 0, 0, 0, 2880*T^3, 10800*T^4, 25920*T^5,  50400*T^6;
               0, 0, 0, 0, 5040*T^4, 20160*T^5, 50400*T^6, 100800*T^7];
        
        Q = blkdiag(Q, Q_k);
        M = blkdiag(M, M_k);
    end
end